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Monte  Carlo  simulations  are  reported  for  1:1  electrolytes  confined  be¬ 
tween  electrically  charged  and  neutral  surfaces.  Both  solvent  and  ions 
are  modeled  as  hard  spheres  with  hard  sphere  and  Coulomb  interac¬ 
tions.  Density  profiles  of  solvent  and  ions  and  electrostatic  potential 
profiles  are  presented  for  electrolyte  concentrations  of  1  and  2M  and 
for  different  cation  sizes  at  300K.  Comparison  of  results  at  charged  and 
neutral  walls  indicates  that  the  density  distribution  of  the  solvent  near 
a  wall  induces  significant  layering  of  ions  and  a  decrease  in  the  potential 
in  the  double  layer.  This  layering  effect  cannot  be  captured  in  simula¬ 
tions  of  the  conventional  primitive  model  (PM)  of  the  double  layer. 


INTRODUCTION 

Systems  with  electrolytes  in  contact  with  a  charged  or  neutral  surfaces  are 
relevant  to  many  areas  of  sciences:  electrochemistry,  biology,  tribology,  and  colloid 
physics  and  chemistry.  Among  the  models  of  the  electrical  double  layer,  the  most 
widely  used  is  the  primitive  model  (PM),  in  which  ions  are  represented  as  point 
charges  that  are  embedded  in  hard  spheres,  the  solvent  is  modeled  as  an  isotropic 
dielectric  continuum,  and  the  surface  is  modeled  as  a  hard  wall  with  uniform  surface 
charge  density.  For  low  surface  charge  and  dilute  electrolyte  solutions,  the  classical 
theory  of  Gouy  (1]  and  Chapman  [2]  (GC)  gives  predictions  of  the  interfacial  ion  and 
potential  distributions  that  are  in  excellent  agreement  with  experiments  and  com¬ 
puter  simulations.  Vigorous  investigation  of  the  PM  by  computer  simulations  [3-9] 
during  the  last  decade  have  revealed  phenomena  of  the  electrical  double  layer  that  are 
not  captured  by  either  the  GC  or  other  theories.  For  high  surface  charge  densities  or 
high  electrolyte  concentrations,  the  ion  density  profiles  in  the  PM  indicate  a  highly 
organized  layering  of  ions  at  the  charged  surface.  For  both  1:1  and  2:2  electrolytes 
at  high  concentrations,  there  is  charge  inversion  with  coion  densities  exceeding  the 
counterion  density  in  the  second  layer.  The  obvious  disadvantage  of  the  PM  is  that 
the  interfacial  structure  of  the  solvent  is  minimized,  although  it  is  well  known  that  ion 
solvation  and  dipole  orientation  at  the  surface,  both  of  which  are  critical  in  determin- 
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ing  the  potential  and  ion  distributions,  are  correlated  with  the  molecular  structure 
of  the  solvent.  The  development  of  more  realistic  models  of  water  and  the  attempt 
to  refine  the  primitive  model  [10]  allow  the  effects  of  the  solvent  structure  on  the 
double  layer  to  be  studied.  The  complexity  of  realistic  models  of  water  limits  studies 
of  these  systems  to  computational  simulations;  however,  very  recent  investigations 
[11-13]  demonstrate  the  importance  of  solvent  structure  on  the  interface  properties. 

In  this  article,  we  present  Monte  Carlo  simulations  of  a  solvent  primitive  model 
(SPM)  which  is  comprised  of  charged  hard  sphere  ions  and  uncharged  hard  sphere 
solvent  molecules.  Ion  density  and  potential  profiles  are  compared  for  electrolytes  in 
contact  with  charged  and  neutral  walls  for  different  electrolyte  concentrations  and 
sizes.  It  is  found  that  the  finite  size  of  the  solvent  results  in  a  highly  ordered  layering 
of  ions  at  both  charged  and  neutral  walls.  For  electrolytes  comprised  of  different  sized 
anions  and  cations,  the  solvent- induced  structure  results  in  a  space  charged  layer  more 
than  two  molecules  thick  at  a  neutral  wall. 


MODEL  SYSTEMS 


The  systems  simulated  in  this  work  are  mixtures  of  1:1  primitive  electrolytes  and 
neutral  hard  spnere  solvent  molecules  contained  between  two  infinite  planar  surfaces. 
One  of  the  surfaces  (at  z  —  0)  is  uniformly  charged  and  the  other  (at  z  =  H)  is 
neutral.  In  the  computations,  the  systems  are  represented  as  a  periodic  array  of 
closed  rectangular  boxes  of  height  H  and  a  square  base  of  length  L  on  a  side.  Let  jV+ 
and  and  q+  and  denote  the  number  and  charge  of  cations  and  anions  in  the 
enclosed  system,  respectively.  Electroneutrality  requires 


1^-9-  -  iV+g+|  =  <xA  , 


(1) 


where  a  is  the  surface  charge  density  and  A  =  L2.  The  interaction  potential  between 
particles  t  and  j  is 


U(rij)  = 


( 


oo  , 

elrt~rd  ’ 


and  between  particle  i  and  the  wall  is, 


rij  <  (di  +  d})/2 
ri}  >  ( di  +  d})/2  , 


(2) 


oo , 

-2  raqiZijt, 


Zi  <  d{/2 
Zi  >  di/2  , 


(3) 


where  Zj  is  the  distance  between  the  wall  and  particle  i.  In  equations  (2)  and  (3), 
di  is  the  diameter  of  particle  i  ,  qi  equals  the  ion  valence  times  the  unit  electronic 
charge  e.  e  is  the  dielectric  constant  and  is  assigned  the  value  78.5  in  our  calculations. 
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SIMULATION 


The  Coulomb  interaction  between  ions  in  a  confined  system  which  is  divided 
into  many  cells  in  the  x  and  y  directions  is, 


U{r)  =  £ 

R 


1 

|R  +  r|  ’ 


(4) 


where 

\R  +  r\  =  [(l  +  x)2  +  (m  +  y)2  +  z2)l'\  (5) 

Using  the  method  initiated  by  van  der  Hoff  and  Benson  [14]  in  the  50’s,  Lekner  [15] 
obtained  the  following  result, 


U  =  4  £  cos(2xlx)Ka\2irl\J(m  4-  y)2  4-  z2}  —  log  [cosh  (2irz)  —  cos(27ry)],  (6) 

lm 


where  K0(x )  is  the  modified  Bessel  function.  For  large  values  of  x,  the  function 
behaves  as  e~x/y/x.  In  the  actual  calculations,  we  only  need  to  retain  terms  corre¬ 
sponding  to  x  <  8,  the  contributions  from  long  range  tail  will  be  of  the  order  ~  10-4. 


In  a  simulation,  we  choose  ionic  concentrations  to  mimick  those  of  the  PM 
electrolyte  in  equilibrium  with  the  bulk  salt  at  a  given  concentration.  First,  we  run 
the  grand  canonical  Monte  Carlo  (GCMC)  simulations  for  the  primitive  electrolyte 
systems  [3-4].  Then  a  canonical  MC  simulation  is  carried  out  for  a  solvent/electrolyte 
system  at  a  total  number  density  of  [5], 

(N+  +N.  +  No)cP/HL 2  ~  0.7,  (7) 

where  N+  and  AL  are  the  average  number  of  free  ions  in  the  confined  PM  fluid.  At 
such  high  particle  densities,  a  grand  canonical  Monte  Carlo  simulation  is  very  hard 
to  perform  due  to  the  difficulties  in  inserting  particles.  This  is  why  a  canonical  MC 
simulation  is  used  in  the  second  stage. 


The  parameters  used  in  our  simulations  are  listed  in  Table  1.  The  probability  of 
accepting  a  move  in  the  SPM  simulation  is  10  percent.  During  the  run,  the  systems 
were  equilibrated  for  10*  steps  and  then  another  (2.5  —  4.0)  x  10®  steps  were  run  to 
accumulate  averages  to  calculate  the  density  profiles  n(x)  and  the  mean  electrostatic 
potential  ip(z)-  Very  long  runs  are  required  compared  to  the  primitive  electrolytes 
simulations,  in  order  to  eliminate  the  noise  in  the  density  profiles.  Taking  the  right 
wall  as  the  reference  point  of  the  potential,  ip(z)  is  obtained  via  the  relation 


Table  1:  Parameters  used  in  this  paper,  cr*  =  ad?  je,  side  length  L,  surface  separation 
H.  Number  of  cations,  anions,  and  neutral  particles  are  N+  .  .V_,  and  V0,  and 
Nc  =  N+  +  N-,  AN  =s  7V_  —  N+  =  a  A,  A  =  L  x  L.  C  is  the  bulk  ion  concentration 
for  PM  electrolytes  in  molar.  The  number  of  steps  run  in  the  simulations  are  given  as 
Nsteps  in  unit  of  106  steps  in  the  last  column.  In  all  the  cases  studied  in  this  paper 
the  neutral  solvent  particle  and  the  anion  have  diameter  of  do  =  d_  =  d  —  4.25  A. 
For  systems  a  ,  6,  and  c,  ln7±  =  -0.127,  0.271,  and  —0.44,  respectively. 


sys 

Hjd 

Ljd 

djd. 

<7* 

Nc 

No 

AN 

C(M) 

Nsteps/ 106 

a 

12.0 

4.36 

1.0 

26 

134 

8 

1.0 

4.0 

b 

12.0 

4.47 

1.0 

32 

136 

14 

1.0 

2.5 

c 

12.0 

4.36 

0.5 

0.42 

26 

134 

8 

1.0 

2.5 

RESULTS  AND  DISCUSSION 

In  the  following  we  report  our  simulation  results  of  the  density  profiles  and  mean 
electrostatic  potentials  for  the  different  systems  that  are  listed  in  table  1.  The  surface 
charge  density,  the  ion  and  solvent  density,  and  the  mean  electrostatic  potential  are 
reported  in  the  dimensionless  forms 

or*  —  Z — (  nT  —  and  ipm  —  (9) 

e  ni0 

where  nlo  is  the  bulk  number  density  for  species  i  and  n,  is  the  number  density  in  SPM. 

To  compare  the  results  obtained  here  with  those  obtained  for  the  primitive  elec¬ 
trolytes,  in  figure  1,  we  show  the  SPM  results  as  continuous,  dotted  or  dashed  curves 
and  the  PM  results  as  open  and  filled  circles  for  the  case  of  cr'  =  0.42  and  C  —  1 M. 

As  is  well  known,  the  Gouy-Chapman  theory  and  simulations  of  the  primitive 
electrolytes  are  in  agreement  for  low  surface  charge  density  and  electrolyte  concen¬ 
tration.  Under  these  conditions,  both  the  density  profile  and  the  mean  electrostatic 
potential  are  monotonic.  In  figure  1,  we  cm  see  that  the  solvent  size  included  in  the 
SPM  simulation  has  a  dramatic  effect  on  the  electrolyte  ion  profiles.  The  presence  of 
solvent  molecules  at  the  walls  induces  strong  structure  in  the  ion  distributions  at  both 
charged  and  neutral  walls.  Near  the  left  wall,  the  neutral  particles  and  counterions 
have  five  layers,  the  coions  show  four  layers  and  near  the  right  wall  there  are  four 
layers  of  each  type  of  particles.  A  small  residue  of  the  coions  near  the  left  wall  is 
seen  in  the  SPM  model  simulation.  These  structural  features  are  absent  in  the  PM 
simulations.  Figure  2  shows  the  mean  electrostatic  potential  for  the  same  system  as 
in  figure  1.  The  potential  decreases  rapidly  near  the  charged  wall  and  decays  to  zero. 
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As  mentioned  above,  the  noise  in  the  density  profiles  prevents  us  from  getting  more 
accurate  results  to  see  a  clear  trend  of  the  potential. 

In  figures  3  and  4,  we  change  the  ion  concentration  to  2 M  and  use  a’  ~  0.4. 
Again  the  layering  of  ions  is  very  strong  and  similar  to  figures  1  and  2.  Furthermore, 
a  charge  inversion  occured  in  the  second  and  third  layer  with  the  coion  density  ex¬ 
ceeding  the  counterion. 

To  examine  the  effects  of  ion  diameters  on  the  highly  structured  ion  distribu¬ 
tions  of  the  SPM,  we  ran  simulations  for  electrolyte  comprised  of  dissimilar  sized 
particles,  namely,  d+  =  0.5d_  .  The  results  are  given  in  figures  5  and  6.  These  results 
show  that  the  coion  peaks  shift  due  to  the  size  effect  but  there  is  still  significant  lay¬ 
ering  structure  in  the  ion  and  solvent  density  distributions.  Ion-size  induced  charge 
separations  are  clearly  s^en  near  the  neutral  wall.  Due  to  the  noise  in  the  density 
profiles,  we  cannot  make  an  accurate  estimate  of  the  potentials  resulting  from  these 
separations. 

From  the  above  analysis,  the  layered  structures  that  appear  in  the  SPM  elec¬ 
trolytes  are  mainly  due  to  the  high  total  fluid  density  resulting  from  the  presence  of 
the  solvent  molecules.  The  small  residue  of  coions  near  the  charged  surfaces  is  due  to 
the  layering  of  the  solvent  and  counterions. 
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Figure  1.  (Upper)  Density  profiles  nm(z)  =  n(z)/no 
Figure  2.  (Lower)  Mean  electrostatic  potential 
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Figure  3.  (Upper)  Density  profiles  n*(z)  =  n(z)/no 
Figure  4.  (Lower)  Mean  electrostatic  potential  i>m{z). 
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Figure  5.  (Upper)  Density  profiles  n*(z)  =  n(z)/n0 
Figure  6.  (Lower)  Mean  electrostatic  potential  0*(z). 


